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1 Introduction 

The spectrum of numerical methods for parabolic evolution equations is extremely 
broad, which attests to the ubiquity and the relevance of such equations. With 
the aim of developing reliable massively parallel algorithms, e.g. for optimiza- 
tion problems constrained by parabolic evolution equations, several attempts have 
been made to go beyond time-stepping methods, see [S', Section 5.1] for a mod- 
est attempt of an overview. In this paper we give a concise Matlab implementa- 
tion of a specific space-time Petrov-Galerkin discretization for parabolic evolution 
equations from 01[3], hoping to provide a basis for possible further developments. 
Spanning just a few lines of Matlab code, it is parallelizable and stable in the 
Petrov-Galerkin sense, which already distinguishes it from conventional methods 
for parabolic evolution equations. Stability implies quasi-optimality of the discrete 
solution in the natural solution spaces, and is an important property in the resolu- 
tion of nonlinear problems. Moreover, the implementation is modular with respect 
to the spatial discretization, admits time-dependent inputs and nonuniform tem- 
poral grids. Since the algorithm is based on an iterative solution of a single linear 
system, another significant advantage to conventional methods is the ability to 
terminate the iteration at a specified global accuracy. 

The model parabolic evolution equation under consideration is presented in 
Section [2] and is restated in a variational formulation. The space-time Petrov- 
Galerkin discrete trial and test spaces that will be used to discretize the variational 
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formulation are introduced in Section [3l In order to obtain stable algorithms we 
develop in Section |3] a generalization of the usual variational framework for lin- 
ear operator equations, called the minimal residual Petrov-Galerkin discretization. 
Choosing bases on the discrete trial and test spaces, it leads to a linear system of 
generalized Gaufi normal equations along with a natural preconditioner. In that 
framework, certain norm-inducing operators play an important role. Specifically 
for the parabolic evolution problem, these are defined in Section [5l In Section 
[6] we detail the Kronecker product structure of the parabolic space-time operator 
and the norm-inducing operators when assembled using space-time tensor product 
bases, and comment on the data structures employed. In the solution process, the 
inverses of the matrix representations (of norm-inducing operators) are required. 
Their Kronecker product structure is discussed in Section The assembly pro- 
cedure for the space-time source is in Section [8l A generalization of the LSQR 
algorithm of Paige and Saunders [16] for the iterative resolution of the linear sys- 
tem is given in Section (9) With those preparations, the Matlab implementation is 
presented in Section llOl and a numerical example is given. We conclude and point 
out some limitations in Section [TTl 

2 Model problem and its space-time variational formulation 

Let D C M", n £ {1,2,3}, be a bounded connected domain with a polyhedral 
boundary F = dD. If n = 1 then D is an open bounded interval; if n = 2 then 
D is a polygon; if n = 3 then _D is a polyhedron. We partition 7^ into two disjoint 
subsets To and Fn such that F = FaYJ Fn- On Fo we will impose homogeneous 
boundary conditions of Dirichlet type. For that reason, the Dirichlet boundary Fq 
is assumed to be of positive measure (with respect to the surface measure which 
will subsequently be denoted by a), i.e., it contains at least one endpoint if n = 1; 
a curve of positive length if n = 2; or a surface of positive surface measure if n = 3. 
Let J = (0,T), T > 0, denote the temporal interval. 

The model for parabolic evolution equations that we consider is the heat equa- 
tion: 

{t,x)eJxD, (1) 
{t,x)eJxFo, (2) 

{t,x)eJxFN, (3) 
{t,x)e{0}xD. (4) 

Here, a, /, g and h are given scalar valued functions, while u is the unknown; 
further, A = ^27=1 'WS^ ^^'^ respectively, denote the Laplace operator on D 
and the derivative in the direction of the outward normal at the Neumann part 
Fn of the boundary. The precise meaning of the heat equation will be fixed by 
means of a well-posed space-time variational formulation in the following. The 
time-independent operator "div(a(a;) gradu(t, a;))" could be replaced by the time- 
dependent one "div(A(t, x)gradu(t,a;))", where A G L°°{J x _D;]R"y^), without 
affecting most considerations below (the technical reason why this is possible is 
given in [131 Lemma 4.4.1]). However, if A is not a finite sum of separable func- 
tions, the implementation becomes significantly less transparent, and we therefore 



dtu{t, x) — div(a(a;) gradu(t, x)) = f{t, x), 

u{t, x) = 0, 

a{x)^it,x) = g(t,x), 
u(t, x) = h(x), 
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discard this case from the onset on; on the other hand, A being a finite sum of 
separable functions entails modifications of secondary relevance to this exposition. 

To motivate the space-time variational formulation of the heat equation, we 
formally test the equation with a function vi on J x D and integrate in space 
and time; the initial condition is tested with V2 on D and integrated in space. 
Integration by parts is performed in space only, and the two resulting conditions 
(one "Vi;i", the other "\/v2") are added together. The solution u will be sought in 
the space X, and the test functions are combined to v = {vi,V2) £ Y := Yi x Y^. 
The spaces X and Y will be specified presently. We write (•,-)-D for the L^{D) 
and the [L^{D)]'^ scalar product, while (•, is the scalar product on L^{rN) for 
the boundary measure a introduced above. Generally, we omit the dependence of 
the integrands on the temporal variable t. In this way we obtain the continuous 
space-time variational formulation 

find ueX such that B{u,v)=b{v) \fv e Y, (5) 

where the system bilinear form B encoding the heat equation is 

B{u,v) := I {dtu,vi)Ddt + j (agradu, gradt;i)ud< + (u(0, •), 1)2)15 (6) 

and the load functional h supplying the source term, as well as boundary and 
initial data, is 

K-") ■■= l^{f,vi)Ddt+ ljg,vi)rj,dt+{h,V2)D. (7) 

Integrating by parts also in time, as in e.g. [6,7J, leads to an alternative space- 
time variational formulation, which, however, will not be discussed below. 
Let us introduce the abbreviations 

V:=Hh,{D), H:=L^{D) and V = [Hh,{D)]' , (8) 

where Hp^{D) denotes the Sobolev space of functions in H^{D) with vanishing 
trace on the Dirichlet boundary Jo, and [Hp^{D)]' its dual. We identify H with 
its dual H' by means of the Riesz isomorphism on H. Then the duality pairing on 
V xV' (or V' X V) is the continuous extension of the if-scalar product onV xV. 
In this way we obtain a so-called Gelfand triple of separable Hilbert spaces 

V ^ H^H' ^V' (9) 

with continuous and dense embeddings. In order for _B:Xxy— >Rtobea 
continuous bilinear form and for 6 : y — > M to be a continuous linear functional, 
it is now natural to take 

X := L^{J;V)r\H^{J;V') and Y := {J ■,V) x H . (10) 

For the definition of the Bochner spaces L^{J\V), H^{J;V'), and the like, we refer 
to e.g. [12] or [151 Chapter 1]. Then dS])-© are well-defined for aU {u,v) e X xY 
whenever 



f eL^{J;[HhAD)]'), geL^iJ-H^/^iFN)) and heL^{D), (11) 
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and 

a£L°°{D) with < essinf a < esssupa < oo. (12) 

For the remainder of the article we assume / e L^{J;L^{D)). The spaces X and 
Y are themselves Banach spaces for the norms ||-||js: and that are given by 

:= ||w||i2(j.y) + ||9tw||i2(j.y,), ueX, (13) 

\\v\\Y:=\\vi\\h^j,y) + \\v2fH, v = {v,,V2)eY. (14) 

We recall from e.g. ,12. Section 5.9.2] or [15, Chapter 1] that any (representant 
of any) u £ X admits a modification on a negligible subset of J such that the 
resulting function coincides with a unique continuous iJ-valued function defined 
on the closed interval J; moreover, the C° norm of the latter is controlled by the 
X norm of the former. In other words, the following embedding is continuous 

X = L'^{J;V)nH\j;V') ^ C°{J;H), (15) 

and in particular 

3O0: \\u{0)\\h <C\\u\\x VweX. (16) 

In this way, the initial value u(0) of any u £ X is well-defined in H. 

With the above assumptions, the space-time variational problem ((5]| has a 
unique solution u & X and the solution depends continuously on the functional 
b G Y' , see [l7i Theorem 5.1] and the references therein. Hence, the solution u also 
depends continuously on the input data /, g and h that define the load functional 
b in ©. 

3 Space-time tensor product discrete trial and test spaces 

The continuous space-time variational formulation ([5]) will be discretized using 
finite-dimensional discrete trial and test spaces Xh C X and Yh C Y built up 
from finite-dimensional "univariate" temporal subspaces E C H^{J), F C L^{J), 
and spatial subspaces Vh C V. These spaces assume the space-time tensor product 
form 

Xh:=E® Vh and Yh := {F ® Vh) x Vh. (17) 

A key feature of this discretization and the implementation given below is the 
modularity with respect to the spatial subspaces Vh ■ 

To specify the temporal subspaces E and F we need to introduce some ter- 
minology. A temporal mesh T is a finite set of points in J = [0,T] containing 
and T. The connected components of J \ T are called the elements of T. Let 
max AT denote the maximal "time-step", i.e., the maximal length of an element 
of T. For a temporal mesh T let T* denote the temporal mesh obtained from T 
by a uniform refinement (each element is split into two smaller elements of equal 
length). Concerning E on the trial side and F on the test side, we will restrict 
ourselves to two types of pairs that differ in the choice of F. 
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Type 1 temporal subspaces. Given a temporal mesh Te, we define E as the 
standard space of continuous piecewise afhne functions, and F as the space of 
piecewise constant functions on Tp := Te- 

Type 2 temporal subspaces. Given a temporal mesh Te, the space E is de- 
fined as above. Let another temporal mesh, 7f be obtained from Te by a 
succession of uniform refinements, i.e., 7f = [T i— > T*]^{Te) for some positive 
n G N. Then F is defined as the space of piecewise constant functions on 7f ■ 

In the first case the discrete variational formulation 

find Uh e Xh such that B{uh,Vh) = biyu) ^vu G Yh (18) 

is an example of continuous Galerkin time-stepping schemes |14| [T] . In the second 
case the dimension of Yh is larger than that of Xh, and the above discrete varia- 
tional formulation is meaningless. A generalization based on residual minimization 
is therefore introduced in Section U) Concerning the stability of the resulting min- 
imal residual Petrov-Galerkin method there is a fundamental difference between 
Type 1 and Type 2 temporal subspaces. This is the subject of the following two 
propositions that summarize the relevant main results from [3l Section 5.2.3]. Note 
carefully that the present concept of stability, namely the validity of the discrete 
inf-sup condition 

. B{uh,vh) ^ „ 

7h := mf sup — f — - — ^ > (19) 

UheXH\{0} y^^YiMO} \\'^^h\\x\\Vh\\Y 

(its role is discussed in Section \^ uniformly in the choice of the temporal dis- 
cretization, is different from e.g. A-stability for time-stepping methods. 

The following measure of self-duality for the spatial subspace Vh C V will be 
needed: 

Kh:= inf sup II y^^'^i'^^'^i • (20) 

Note that Kh is bounded (independently of Vh) and necessarily positive for a 
finite-dimensional Vh- 

Proposition 1 Let {0} ^ Vh G V be a finite- dimensional subspace. Let E C 
H (J) and F (Z L (J) be of Type 2. Then there exists a constant 7q > indepen- 
dent of Vh, E and F, such that the discrete inf-sup condition (|19p holds for the 
discrete trial and test spaces (|17p with 

Ih > jQi^h- (21) 

We remark that 70 in (|2ip , as a function of the number of refinements between 
Te and 7f, is monotonically increasing and saturates (exponentially quickly). 

Type 1 temporal subspaces, on the other hand, do not lead to unconditional 
stability of the form (|2ip . To formalize this, we define the CFL number 

CFU:=maxZ\rE sup JMt^- (22) 
Xhev,A{o} WXhWv 
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Proposition 2 Let {0} Vh C V be a finite- dimensional subspace. Let E C 
H^{J) and F C L^{J) be of Type 1. Then there exists a constant 70 > indepen- 
dent of Vh, E and F, such that the discrete inf-sup condition p9p holds for the 
discrete trial and test spaces (|17p with 

7/^>70K/^min{l,CFL-l}. (23) 

In general, the dependence on the CFL number cannot be improved. 



4 Minimal residual Petrov-Galerkin discretization 

In this section we consider an abstract continuous bilinear form B : X xY 
where X and Y are Hilbert spaces with norms and \\-\\y ■ Let 

iiDii B(w,v) 
\\B\\ := sup sup 



vex\{o} v^Y\{o} \\w\\x\\v\\y 

denote the norm of B. Further, let 6 be a linear continuous functional of Y . For the 
remainder of the section, two finite-dimensional subspaces X^ C X and Y^ C Y 
are fixed. We aim at relaxing the discrete variational formulation (|18|) to admit the 
case dimX^ < dim Y/i. To guarantee well-posedness, the discrete inf-sup condition 
of B on Xh X Yh will be essential (cf. Proposition [1]): 

Ih ■■= mf sup J. — \. — — {r— > 0. (24) 

For practical reasons we introduce norms |||-|||x and |||-|||k defined on X^ and Y^ 
that are induced by (positive definite) linear continuous operators M : X^ — > X' 
and N : Yh ^ Y' as follows: 

:= iMwh)iwh), Wh e Xh, (25) 
llvhlfy ■■= iNvh){vh), VheYh. (26) 

The operators M and A'' are moreover assumed to be symmetric, i.e., {Muh){wh) = 
{Mwh){uh) for all Wh,Uh G X^, and similarly for A^. Let < c^m < Dm and 
< d^v < L)if be constants such that 

cImHx <\\-\\x <DmI\Ix and d^vlMlk < IMIv < ^7v||M||r (27) 

on Xh and Yh, respectively. We emphasize that the operators M and A'', the 
induced norms, and hence the constants in (|27p may depend on h, but in this 
section, Xh x Yh is being considered fixed to lighten the notation. 

Instead of the usual discrete variational formulation we now introduce the 
discrete (functional) residual minimization problem 

Uh ■■= argminRhiwh), (28) 



where 



D f ^ \B{wh,Vh) - b{vh)\ ^ ^ 

Rh{Wh) ■■= sup ^ '-, WheXh, (29) 

vi,eY^\{0} \\\'"h\\\Y 



is the (functional) residual. The following can be shown [3] . 
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Theorem 1 Let the discrete inf-sup condition (|24p hold. Then there exists a 
unique Uh £ Xh for which 

Rh{uh) < Rh{wh) e Xh (30) 

holds. Moreover, satisfies the quasi- optimality estimate 

\\u-Uh\\x<Ch inf \\u-wj,\\x with Ch=^^^^ (31) 
for any solution u £ X of "B(u, ■) = b on Y". 

Proof (Sketch) Invoking the open mapping theorem, one can show that the the 
mapping 6 i— > is weU-defined, hnear and continuous, and its norm is dominated 
by Thus the composition u n- Bu i— > it^ is a continuous projection with 

norm not exceeding Ch given in (|3ip . An application of [19, Lemma 5] finishes the 
argument. 

Having introduced the somewhat nonstandard variational definition (I28p of the 
discrete solution, let us describe its computable algebraic equivalent. To that end 
let ^ C Xh and tf' C be bases for the respective discrete spaces. Having fixed 
the pair Xh x Y^, the possible dependence on h is again omitted. The algebraic 
representants of the system bilinear form B, of the load functional b and of the 
norm- inducing operators A'' and M, are defined with respect to the chosen basis 
in the usual way, 

B:=B(<P,tf'), b:=6(if), N := (Artf')(ip'), M := (Af -?)(<?), (32) 

or in componentwise notation B^^ = B{(j),ilj), = 6(V'), N^'^ = {Nil;){ip'), 
= (M(^)((^') for e <l> and ^ The matrix B is injective if and 

only if the discrete inf-sup condition (|24p holds; further, N and M are symmetric 
positive definite matrices due to the analogous properties (|27p of the operators N 
and M. Using these definitions, the discrete functional residual minimization (|28p 
can be seen to be equivalent to the discrete algebraic residual minimization 

u := argmin ||Bw — b||N-i. (33) 

weR* 

A vector u is a stationary point of (|33p if and only if it satisfies the first order 
conditions, namely the Gaufi normal equations 

b"^N"^Bu = B"^N"^b. (34) 

Since we residual error is measured in the ||-||n-i norm, we will refer to p4p as 
the generalized Gaufi normal equations. 

If the discrete inf-sup condition (|24p holds, the matrix B^N~^B is symmetric 
positive definite; then, the Gaufi normal equations (|34p . and therefore also the 
discrete minimization problem (|33p . have a unique solution. Finally, the matrices 
M and B^N~^B are spectrally equivalent with Section 4.1] 

lhdMdN\\M\M. < ||w||btn-ib < ||B||i^A/^Jv||w||M Vw e M*. (35) 

Therefore, M is a natural preconditioner for the Gaufi normal equations (|34|1 . By 
(|35p . the quality of this preconditioner is controlled by the discrete inf-sup constant 
Ih in (|24p and the norm equivalence constants in (|27p . and does not depend on 
the choice of the basis. 



8 



Roman Andreev 



5 Parabolic space-time preconditioners 

In Section |4] we admitted general norm-inducing operators M and A'' on the fixed 
pair of finite-dimensional discrete trial and test spaces Xh x Yh- For the space- 
time variational formulation ([Sj several practical choice are available [3] Chapter 
6]. Here, to simplify the exposition, we will only use the canonical choice of the 
Riesz mappings defined (on all of X and Y) by 

{Mw){w):=\\w\\l = \\w\\l.^j.y^ + \\dtw\\hf^j.y,), wEX, (36) 

and 

{Nv)iv):=\\v\\^ = \\vi\\h^j.y) + \\v2fH, v = {vi,V2)eY. (37) 

These definitions extend to the off-diagonal by the imposed symmetry of M and 
A'', but this will not play a role in the following. In these formulas, the V = Hp^ {D) 
and the V' = [Hp^iD)]' norms are taken to be the "energy norms": 

\\(j\\v ■■= (a grad (T, grader) D, a e V, (38) 
yfv,:={^,A-\ip)D, ipeH = L\D), (39) 

where A is the operator A : V ^ V' , a i— > (agradu, grad-)^!, and up G V' is 
the functional on V defined by Lip := {ip, ■)d, f £ H. For a functional /3 G V' , 
the definition extends to ||/3||v' := f3{A~^l3), but this will be irrelevant in the 
implementation below. 



6 Kronecker product structure of the discretized operators 

6.1 The system bilinear form 

Recall from Section [2] the definition of the system bilinear form 

B{u,v) := J {dtu,vi)Ddt + J {a giadu, giadvi) odt + {u{0,-),V2)d 

for the space-time variational formulation of the model parabolic evolution equa- 
tion, where ueX = (J;V) n H\J;V') and v = {vi,V2) eY = L^{J\V) x H. 
As described in Section [S] we consider two types of discrete trial and test spaces. 
In either case these have the form 

Xh=E®Vh(lX and Yh = [F ®Vh) y. Vu dY, (40) 

where E C H^{J) is the space of continuous piecewise affine functions on a tem- 
poral mesh, F C L^{J) is the space of piecewise constant functions on the same 
mesh (Type 1) or on its 7i-fold uniform refinement (Type 2), and Vu <Z V \s & 
finite-dimensional subspace. 

For the remainder of the section we fix the spatial discretization Vh C V with 
a basis S C Vh, and the temporal mesh Te = {0 = to < ti < . . . < Ik = T}. Let 
7f = {0 = t'o < t'l < . . . < t'l^i = T} be either Te or any n-fold uniform refinement 
of Te ■ Let E be as above, and let F denote the space of piecewise constant functions 
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with respect to the temporal mesh 7f ^ Te , which possibly refines Te • As basis 
for E we take the usual hat functions := {9k : k = 0, . . . , K} C E defined by 
^fe(^fc) ~ ^fcfc- particular, the only function that does not vanish at t = is ^?o- As 
basis for F we take the indicator functions S := {^^ := X(t'j. ^.tj.) : fc = 1, • • • , K'} 
on the elements (ifc_i,tfc) of the temporal mesh 7f ^ Te- These univariate bases 
are first combined to the collections <P C X , C Yi and ^2 C Y2 as 

4>:={e(x)a:eee,aeIJ}, 'Pi := {^®a : 9 e S,a e U}. (41) 



These now yield bases 



<PcXh and 1^ := {<Fi x {0}) U ({0} x S) C Yh (42) 

for the discrete trial and test spaces and Yh. Discretizing the bilinear form B 
using these tensor product bases <P and ^ as described in abstract terms in Section 
m leads to 



ef ® Mx 



(43) 



where a) the "temporal FEM" matrices Cf^,Mf-^ e R^'^® and the row vector 
ef G , have the components 

[CrUe = j^9'{t)at)dt, [M[%e = j mmdt, [ef]e = Sg^g, (44) 

with the prime denoting the derivative with respect to t, and b) the usual "spatial 
FEM" mass and stiffness matrices M^; , A^ £ R"^ ^ ^ are given by 

[M^^J^cr = / a{x)a{x)dx, [Ax\aa = / a(a;) grad(T(a::) • grada(a:)da;. (45) 

J D J D 

Let us comment on the assembly of the temporal FEM matrices. First, if 
Te=Tf (Type 1) then [C[%g e {1,-1,0} and [M[%g G {i|/|,0} depending 
on whether ^ and 9 are both nonzero on the temporal element / of Te having length 
|/|, and on the sign of 9' there. Therefore, assume now that Tf is obtained from Te 
by a succession of uniform refinements (Type 2). Let Te denote the first uniform 
refinement of Te, and let E* be the space of continuous piecewise affine function 
on Te with the hat function basis O* . Let Ct^ and ^ denote the matrices 
as above, but with Te in place of Te- Consider now the embedding operator 

: E ^ E*. Define the prolongation matrix [S^]g*g by S^9 = Ee-ee* Sf.961*, 
9eO. Then 

Cf^ = Cf^'sf and = Mf^*Sf . (46) 

Moreover, denoting by G T the node for which 9(tg) = 1, and similarly for 
tg* G r*, we have "Te = S^Te" , i.e., 

tg. = J2[^^]o'etg y9*ee*. (47) 

eee 
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6.2 The norm-inducing operators 

With the norms on V and V' taken to be (|38p -(|39 p . the discretized operators M 
and A'^ defined in Section [5] assume the form 

M = Mf (»A, + Af ®(M,A-1m^) (48) 

with the "temporal mass and stifl^ness" matrices 

i^fUe = 1^ Omm, [Af = ^ 9'{t)e'{t)dt, (49) 



and 



Mf ® A^ 



J 



with the "temporal mass matrix" 

[Mf = / mmdt- (51) 

6.3 Data structures 

The Kronecker product structure of these matrices suggests regarding a vector 
w e R^^® as a rectangular array with rows and #6* columns. Let Vec(w) 
denote the "vectorization" of such an array, i.e., its columns are collected one after 
another into one column vector Vec(w) of length #(£■ x 0). Now, if T £ R®^® 
and X e R^^'^ are matrices then 

(T ® X) Vec(w) = Vec(XwT"^). (52) 

In the implementation we will exclusively use the representation as rectangular 
arrays. Moreover, load vectors derived from load functionals d £ Y' will be stored 
as pairs d = (d^, d^) with d^ £ R'^^" (rectangular array with #17 rows and 
columns) and d^ £ R^ (column vector of length i^S) in the form of a Matlab 
structure {dl, d2}. To these, a formula analogous to (|52p applies. In particular, we 
never store the operators B, M and N (or its inverses) as matrices. 

7 Inverses of the space-time parabolic preconditioners 

Consider the matrix representations M and N of the space-time parabolic pre- 
conditioners given in (|48p and (|50|1 in Section [6l In order to solve the generalized 
Gaufi normal equations (|34|) with M as preconditioner, we need to (approximately) 
compute the inverses and N~^. 
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7.1 The test side 

For N we simply use the representation 

which again has Kronecker product structure. For problems with a large number 
of spatial degrees of freedom, the inverse A^^ may replaced by any approximate 
inverse, e.g. several cycles of a multigrid method. 




7.2 The trial side 

The representation of requires more discussion, as it is not (obviously) of 

Kronecker product structure. We will obtain a simplified expression for M"^ by 
diagonalizing Nlf and Af'. For this discussion, let us drop the superscript (•)^. 
Consider therefore the generalized eigenvalue problem of finding v G and A G R 
such that Atv = AMtv (in place of Mt, one could use the mass-lumped version of 
Mt , or simply the diagonal matrix that coincides with Mt on the diagonal) . Let It 
denote the identity matrix of the same size as Mt and At. Since At is symmetric 
positive semi-definite and Mt is symmetric positive definite, all eigenvalues are 
nonnegative and the eigenvectors may be chosen to form an Mt-orthogonal basis: 
there exists a square matrix Vt collecting the eigenvectors in its columns, and a 
diagonal matrix Df containing the eigenvalues on the diagonal, such that 

VjMtVt = It and AtVt = MtVtDt. (54) 

Let us set Tt := MtVt. The first identity implies V^"^ = VjMt = TJ , which 
may be used to verify 

Mt = TtTj and At = TtDtXj. (55) 
Inserting these into the expression (|48p for M we find 

M = (Tt ® Ix)(It O A, Dt O (M,Aj'M,))(Tr Ix) (56) 
and therefore 

M-i = (Vt Ix)(It ® A, -h Dt {M^A-^M^))-\VJ ® I,). (57) 

Let je be the square root of the entry of Dt on the diagonal in position 6, i.e., 
Dt = diag((7|)eGe)- Now, recall from Section|6]the convention that w G R^^® is 
stored as a rectangular array with #0 columns. Applying M~^ to such a vector 
w will be done as follows 

1. Compute w^^' :=wVt. 

2. For each column w^^' of w'-^' compute the corresponding column Wg^' of w^^' 

by 

w'"' := (A, +7|(M,A-^M,))-V^^\ (58) 



3. Compute w^^^ := w'-^WJ . 
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Then M = w^'^'. The computation (|58p can be done in parallel over the 
columns; we will further make use of the identity 

(A, + 7i(M,A-iM,))-i = H;1A,H^i = Hi^A.H^i (59) 

where H± := A^; ± i'ygls/lx is precisely the FEM discretization of the Helmholtz 
operator [A ± Id) with imaginary frequency ±17©. In fact, since the input is 
real, it would suffice to use ReHjj|^ instead of (|59p . Interestingly, such Helmholtz 
problems appeared in the context of (parabolic) evolution equations in [18l[8] , but 
in the present case only in the representation of the preconditioner M~^. These 
may therefore be inverted approximately. 

8 Assembly of the space-time load vector 

As in the previous section, let F be the space of piecewise constant functions on 
a temporal mesh Tf ; let S be the basis on F consisting of indicator functions on 
the elements of Tf; finally, 17 C is a basis for a finite-dimensional subspaces 
Vh C V = Hp^{D). The load functional b{vi,V2) defined in Section [2] can be 
rewritten as 

b{vi,0) + b{0,V2) = 1^ {(/, vi)d + (5, vi)rjdt + {h,V2)D, (60) 

for {vi,V2) e Y = L^{J\V) X H, where h e H = L^{D), f G L^{J;H) and 
g G L'^{J\H^^'^{rN)) are given functions. Accordingly, the load vector b that we 
obtain as outlined in Section|3]consists of two parts, say, b^ G R'^^" and b^ G R'^. 
The second part is given in componentwise notation by 

[b']^ = (ft, a)D + (0, (T)r« , aeS, (61) 

which is the usual spatial FEM load vector for the function h G L^{D) with 
zero Neumann data. In the remainder of this section we describe how b^ may be 
obtained from the usual spatial FEM load vectors. Observe first that, fixing an 
indicator function ^ G supported on the closed interval /, we may employ a 
quadrature rule on / to define 

[h\^ ■.= Y.wi[{f{ti,-),a)D + {g(.ti,-),a)r^] (62) 

~ l^{{f,vi)D + ig,vi)r^}dt, aeS, (63) 

where & I are the quadrature nodes and G M are the quadrature weights 
(equal to zero for r G N large enough). Now, each term in the curly brackets {. . .} 
is a load vector for the function f{ti,-) with Neumann data g{tl,-), which can 
be assembled using standard spatial FEM routines. This presupposes sufficient 
regularity of the functions t 1— > f{t, •) and t 1— g(t, •) on the interval /. The com- 
putation of individual columns of b (i.e., for each given ^ & S) can be performed 
in parallel. 

If Tf = Te then using the trapezoidal rule on each temporal element / leads 
to a system Bu = b that admits a unique solution, which on the nodes of Te 
coincides with the solution obtained by the Crank-Nicolson time-stepping method 
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9 Generalized LSQR algorithm 

In Section |4] the minimal residual Petrov-Galerkin discretization was shown to 
lead to a system of generalized Gaufi normal equations. Efficient iterative Krylov 
subspace methods for such problems are known One option is the LSQR 

algorithm due to Paige and Saunders [16] applied directly to the preconditioned 
equation B"^BQ = B"^b with B = N'^^'^BM.-^/'^ , u = M'^/^u, b = N'^/^b, 
where M"^''^ and N~^/^ denote the inverses of the square roots of the (s.p.d.) 
matrices M and N. We now reformulate this version it in such a way that only 
the inverses and need to be applied, but not the square roots, cf. [9]. 

To compute an approximate solution Ui* w u to the generalized Gaufi normal 
equations (|34p using M as a preconditioner, the algorithm proceeds as follows: 

1. Initialize 

(a) do := 

(b) (vo,vo,/3o) := NORMALIZE(b,N) 

(c) (wo,wo,Qo) := Normalize (B"''vo, M) 

(d) po := ||(«o,/3o)||2 

(e) uo := 

(f) 5i = ao, 71 = 

2. For i = 1,2, . . . ,i* do the following steps (until convergence) 

(a) di := Wi_i - (ai_i/3i_i/pf_i)di_i 

(b) (vi,Vi,/3i) := NORMALiZE(Bwi_i -aj_iVi_i,N) 

(c) (wi,w,,Qi) := Normalize(B''"vi - /3iWi_i, M) 

(d) :=||(5.,/3.)||2, 

(e) Ui := u,_i + ((5i7i/pi )d, 

(f) 5i+i := -5iai/pi, -ji+i := ^iPi/pi 

Here, Normalize : (s,S) ^ {z,z,z), with S s.p.d., is the procedure: 

1. Solve Ss = s for s 

2. Set z := x/s^s and (z, z) := {z~^s, z~^s) 

As long as the order of the statements is unchanged, the subscripts etc., 
can be ignored in the implementation. In our implementation we will limit the 
number of iterations and allow the iteration to exit when the normal equations 
residual ||B"^Bu - B'^b||2 = ||B'^N"^Bui - B'^N"^b||M-i falls below a tolerance 
threshold. This residual is available for each i = 0, 1, . . . as 7^-1-1 following 

step (f). See [10] for further discussion of stopping criteria for the LSQR algorithm. 



10 Overview of the Matlab code 

In the code given below the naming convention parallels that of Sectionsl6}]8] Thus, 
the mesh Tf is called tf, the matrix Mf^ is called MtFE, the vector u is called 
u, and so on. The subroutines that are related to the temporal FEM are prefixed 
with "femT_", those related to spatial FEM with "femX_". The only subroutine 
that mixes temporal and spatial FEM is femTX.assemLoad for the assembly of the 
space-time load vector b. 
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10.1 Main file 

We provide the commented code for the main file. The code is embedded in a 
Matlab function spacetime. 

1 function spacetime 

Initialize the spatial FEM (load the mesh, etc., here into global variables) and 
compute the spatial FEM mass and stiffness matrices. The flag true indicates that 
this is the first-time initialization. 

2 f emX_init (true ) 

3 [Mx , Ax] = f emX_MA () ; 

Define a temporal mesh with K elements on the interval J = (0, T) with T = A: 

4 T = 4; K = 100; 

5 TE = T * sort([0, rand(l, K-1), 1]); 

Determine the number of uniform refinements to go from Te to Tf, see Section |6l 
Setting use_mrpg = true amounts to one refinement. 

6 use_inrpg = true ; 

Compute the refined mesh Tf and the temporal FEM matrices as described in 
Section [6} 

7 [MtFE , CtFE , TF] = f emT_assemFE (TE , use.mrpg); 

8 [MtE , AtE] = f eniT_asseniE (TE) ; 

9 MtF = f eniT_asseniF (TF) ; 

Define the function that computes w n- Bw using the matrix representation given 
in Section [6] 

10 function Bw = B(w) 

11 Bw = { Mx * w * CtFE' + Ax * w * MtFE', Mx * w(:,l) }; 

12 end 

Define the function that computes v n- B^v using the matrix representation given 
in Section [6l 

13 function Btv = Bt(v) 

14 Btv = Mx ' * v{l} * CtFE + Ax' * v{l} * MtFE; 

15 Btv(:,l) = Btv(:,l) + Mx ' * v{2}; 

16 end 

Define the function that computes d i— > N~^d with as in Section[71 

17 function iNd = iN(d) 

18 iNd = { Ax \ (d{l} / MtF), Mx \ d{2} }; 

19 end 

Define the function that computes w n- M~^w using the algorithm and formula 
(|59p given in Section[71 Two comments are in order. First, the symmetric positive 
semi-definite matrix A.f is singular, possibly leading to a small negative approx- 
imately computed eigenvalue. Before taking the square root we therefore round 
negative eigenvalues to zero. Second, the result of an application of (|59|1 to a real 
vector is again a real vector, which is enforced by taking the real part. This trun- 
cates the round-off error accumulated in the imaginary part and reestablishes the 
data type of reals. The loop can be performed in parallel. 
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20 
21 
22 
23 
24 
25 
26 
27 
28 
29 



[VtE , DtE] = eig (f ull (AtE) , full(MtE)); 
gamma = sqrt (max (0 , di ag ( DtE ) ) ) ; 
function iMw = iM(w) 



iMw = w * VtE ; 



parfor j = 1: length (TE) 



H = Ax + li * gamma(j) * Mx ; 



iMw(:,j) = real(H \ conj(Ax * (H \ iMw ( : , j ) ) ) ) ; 



end 



iMw = iMw * VtE ' ; 
end 



Compute the space-time load vector b using a quadrature rule according to Section 
|8l The Matlab functions f , g and h are assumed to be available, e.g. as Matlab 
files in the same directory. Further, QR_Trapz is the trapezoidal quadrature rule as 
described in Section [1021 

30 b = f emTX_assemLoad(TF , Of, Og , Oh, QR.TrapzO); 

Set the tolerance and the maximal number of iterations for the generalized LSQR 
algorithm from Section[9]and run it on the Gaufi normal equations (|34[) with M as 
preconditioner. The solver may provide additional diagnostic output parameters 
that are ignored here. 

31 tol = le-8; maxit = 100; 

32 u = glsqr((SB, SBt , b, tol, maxit, aiM , (SiN); 

Finally, plot several temporal snapshots (equispaced in time) of the numerical 
solution. These are obtained by linear interpolation of the values at time nodes 
and stored in the array U. 

33 U = interpl(TE, u', linspace (min (TE) , max(TE), 5))'; 

34 for k = 1 : size (U , 2) 

35 subplotd, size(U,2), k) ; f emX_show (U ( : , k) ) ; 

36 end 

Here the Matlab function spacetime ends. 

37 end 

10.2 Assembly of the space-time load vector 

The assembly of the space-time load vector b is performed in the Matlab function 
femTX_assemLoad. It receives the mesh 7f as a (row) vector tf, as well as function 
handles f, g and h. The function handles f and g, when called with one argu- 
ment, say tr, are expected to return function handles to functions that depend on 
the spatial variable only and describe f{tr, •) and g{tr, •). The function handle h 
describes the initial condition h{-). Finally, QuadRuie is a function handle that re- 
ceives a 2-component vector describing the endpoints of an interval / and returns 
a quadrature rule on / in the form of a vector of nodes {ti)r=i,...,R and a vector 
of corresponding weights (wr)r=i,...,H, as well as the number of nodes R. 

1 function b = f emTX_ as s emLo ad ( TF , f , g, h, QuadRule ) 

Assemble the part b^ of the load vector from the initial condition h{-). 

2 b2 = femX_b(h, ( varargin ) 0) ; 
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Assemble b by iterating over the temporal elements defined by Tf as described 
in Section [51 The outer loop can be computed in parallel. 

3 bl = zeros (length (b2) , length ( TF ) -1) ; 

4 parfor k = l:size(bl,2) 

5 [tl , wl , R] = QuadRule (TF ( [k k + 1])); 

6 for r = 1:R 

7 bl(:,k) = bl(:,k) + wl(r) * f emX.b (f (tl(r ) ) , g(tl(r))); 

8 end 

9 end 

The assembled vectors and b^ are combined into a Matlab structure. 

10 b = {bl , b2}; 

Here the Matlab function femTX.assemLoad ends. 

11 end 



10.3 Assembly of the temporal FEM matrices 1 

Let us comment in some detail on the computation of the temporal FEM matrices 
and Mf ^ by means of the Matlab function f emT.assemFE. It receives a tempo- 
ral mesh To ^ Te and the number nref of uniform refinements to be performed on 
To to obtain Tf- If nref is interpreted as one or zero if it has the logical value true 
or false, respectively. The output is Cf ^ and Mf^, and the nref -fold refinement 
of To stored again in the variable TO. 

1 function [MtFE , CtFE , TO] = f emT_assemFE (TO , nref) 

If no refinement is to be performed, the temporal FEM matrices Cf ^ and Mf^ 
can be computed directly. We remark that temporal meshes are stored as row 
vectors. 

2 K = length (TO) ; 

3 if (nref == 0) 

4 MtFE = spdiags (dif f (TO) ' * [1/2 1/2], 0:1, K-1, K); 

5 CtFE = dif f (speye (K) ) ; 

6 return 

7 end 

Otherwise, the prolongation matrix Sf" is computed first, see Section |6. II As can 
be seen from (|47p . it coincides with the matrix representation of an interpolation 
operator (a more efficient but lengthier implementation is possible here) . 

8 StE = sparse (interpl (1 : K , eye(K), 1:(1/2):K)); 

Perform a uniform refinement of the current mesh To and pass the new mesh 
recursively to f emT_assemFE. Hence, the number of refinements still to be performed 
decreases by one. We obtain the matrices ^ and Cf^ with respect to the 
meshes Tf and To*, and the nref-th refinement of To (which is the desired Tf). 

9 [MtFEs , CtFEs , TO] = f emT.assemFE ( ( StE * TO')', nref-1); 

Apply the prolongation matrix according to (|46p to obtain Mf^ and Cf^. 

10 MtFE = MtFEs * StE; CtFE = CtFEs * StE; 

Here the Matlab function femT.assemFE ends. 



11 end 
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10.4 Assembly of the temporal FEM matrices 2 

The computation of the temporal stiffness and mass matrices Af, Mf' and is 
a routine task, and the Matlab code is therefore provided for completeness without 
further comment. 

1 function [MtE , AtE] = f eiiiT_assemE (TE) 

2 K = length(TE); h = diff(TE); g = l./h; 

3 MtE = spdiags([h 0; h] ' * [1 2 0; 2 1] /6 , -1:1, K, K); 

4 AtE = spdiags([g 0; g] ' * [-1 1 0; 1 -1], -1:1, K, K); 

5 end 

1 function MtF = f eniT_asseniF (TF) 

2 K = length (TF) -1; 

3 MtF = sparse(l:K, 1:K, abs ( dif f ( TF ) ) ) ; 

4 end 



10.5 Example 

Let us set a = I for the heat conduction coefHcient and f(t,x) = sm(t) for the 
source term in ([1]), g = for the Neumann data (|3]) and h = for the initial 
condition Using the well-documented open-box 2d spatial FEM discretiza- 
tion from '2', Sections 2-8] the code described in this section (and some graphical 
postprocessing) produces a figure similar to Fig. [T] 




Fig. 1 Snapshots of the solution at t = 0, 1, 2, 3, 5, produced by the code given in Section [TOl 
Dark corresponds to low, bright to high values. 



Continuing with these parameters, we set T = 20. In Fig.[2]we document the ac- 
curacy of the discrete solution and the number of iterations of the generalized least 
squares solver for equidistant temporal meshes Te of different size and different 
numbers of refinements between Te and Tf • We observe that for Type 1 temporal 
subspaces (no refinement; this method is equivalent to the Crank-Nicolson time- 
stepping), the number of iterations first increases with #Te but then decreases 
again. For Type 2 temporal subspaces (one or more refinements) , the number of it- 
erations varies with #Te much less and is consistently smaller. Taking the number 
of iterations as an indicator of the conditioning of the generalized normal equa- 
tions with M as preconditioner and taking (|35|) into account, these observations 
are in agreement with Propositions [1] and [2] Replacing the trapezoidal rule by a 
higher-order quadrature does not significantly change that picture. Note that this 
is not the number of iterations need to reach the accuracy of the same order of 
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magnitude as admitted by the temporal discretization, but the number of itera- 
tions needed to achieve a fixed small residual in the preconditioned Gaufi normal 
equations as described in Section |21 
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Fig. 2 Accuracy of the discrete solution in the space-time norm (left) and the number 

of generalized LSQR iterations (right) as a function of the number of temporal elements, see 
Section nTJ31 



11 Conclusions 

A concise Matlab implementation of a space-time simultaneous discretization and 
solution algorithm for parabolic evolution equations that is stable, parallelizable, 
and modular has been presented. It admits nonuniform temporal meshes and time- 
dependent inputs. Very efficient natural preconditioners for the iterative resolution 
of the resulting single linear system of equations are available. Extensions to higher 
order in time and to space-time compressive algorithms are possible 4,3,5 . 

Let us point out what we believe is currently the major obstruction to massively 
parallel computations based on the algorithm presented here. Allowing arbitrary 
temporal meshes, the simultaneous diagonalization of the temporal FEM matri- 
ces discussed in Section [3 leads to a nonlocal-in-time transformation Vt . If the 
solution vector is split over multiple computational nodes along the temporal di- 
mension, which is natural, the application of this transformation may result in 
heavy communication. If, however, the temporal mesh is uniform, or derives from 
a successive dyadic partition of a coarse mesh, it may be possible to use more 
efficient fast Fourier or wavelet based transforms for this particular step instead. 
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